This workbook centralises analyses of the PHQ-9, GAD-7 and MASQ total scores longitudinal analyses in the RAMP and COPING samples - identifying the best fit average group trajectories (Latent class growth analyses; LCGA). This corresponds with analytic step 1 in the pre-registration for project “Moderators of symptom trajectories of depression and anxiety during the COVID-19 pandemic”.
The analytic steps that will be followed for PHQ, GAD and MASQ are:
Step 1. Identifying overall group trajectories of anxiety and depression (Latent Growth Curve Modelling; using total scores across all time points for whole sample) Step 2. Assess fit and stability of models Step 3. Identify and conclude the best fitting model Step 4. visualise all relevant trajectories alongside observed data-based trajectories.
Setup
note: some of the install and update processes require command line verification in the console. If you are running this script for the first time and it does not proceed, open and check the console
colour paletttes and plot label list
RAMP colours (first) from website and second is a generated left tetradic (rectangular) colour complement of the RAMP logo green for contrast and complement.
Identify the best fitting latent growth curve form, i.e. that which fits all of the data best, on average, for each outcome.
The MPlus input time interval specifications go from increasing by 1 for the first 8 follow-ups (we measured with two week intervals) to increasing by 2 for the remaining timepoints when we moved to monthly. The choice of interval is arbitrary, but the relative difference should be retained.
We dropped timepoints 1 & 2 (see ‘Overview of model alterations and adaptations’ below)
Four pieces
These dates are based on guidelines in England - TO CHECK no major discrepancies across Scotland; NI; Wales
Slope includes data points: Baseline (7 Apr 2020) - Follow-up 5 (16 Jun 2020)
Slope includes data points: Follow-up 6 (30 Jun 2020) - Follow-up 11 (20 Oct 2020)
Slope includes data points: Follow-up 12 (17 Nov 2020) - Follow-up 15 (9 Feb 2021)
Slope includes data points: Follow-up 16 (9 Mar 2021) - Follow-up 17 (6 Apr 2021)
Notes on lockdowns
England:
Lockdown 2.0: 31 Oct 2020-1 Dec 2020
Easing of restrictions + tier system: 2 Dec 2020 (Tier 3) - 20 Dec 2020 (Tier 4 introduced)
Lockdown 3: 6 Jan 2021 - 8 Mar 2021
Scotland:
Wales:
Northern Ireland:
After manual checks of each stage, an issue with covariance coverage was detected. The resolution is detailed here.
The covariance coverage (% of participants common to any two time points) between some time points is lower than MPlus default (0.1; 10%). This prevents algorithm for missing data being applied (MLR; maximum likelihood with robust standard errors) and standard errors cannot be computed, even when setting the coverage expectations to zero.
Looking at our data structure and MPlus output, this is due to missing data for follow-up 2 & 3 (RAMP only, no COPING data) resulting in perfect zero covariance coverage with follow-up 17 (currently COPING only; RAMP data extraction issue).
CHANGE TO MODEL: Now running without timepoints 2 & 3
phq_check <- readRDS("../../../data_clean/phq/phq.clean_merged_total_scores.rds")
# get a list of the total score columns
total_cols <- grep("total_Wave", names(phq_check))
#count how many columns we have
n_cols <- length(total_cols)
# create a list of how many NAs are missing across these 18 columns per row
list_n_missing <- apply(phq_check[total_cols],1,function(x) sum(is.na(x)))
# count how many people are missing data for all columns (NA = total_cols)
list_all_missing <- sum(list_n_missing ==n_cols)
print(dim(phq_check)[1])[1] 52416
[1] 10410
Important note: all numbers correspond exactly with MPlus output numbers (n per wave and n missing all observations etc).
Don’t rerun models if there is already output for it (to save computational time), do show the output in the console IN FINAL MAKE SURE SCRIPT STATES eval=FALSE so that MPlus models aren’t rerun
This reads in all output, including the gh5 plot files
at the moment this doesnt say which model the warning is from, so any warning requires checking to identify its source. Though order in list will give some indication
[1] “c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10410", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10410", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("All continuous latent variable covariances involving LIN4 have been fixed to 0", "because the variance of LIN4 is fixed at 0.") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10410", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10410", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10410", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS")”
[2] “c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10410", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10410", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10410", "2 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10410", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10410", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10410", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS")”
Warnings are all about missing data (missing on all obs, this number lines up with our expectations from the dataset in the section above)
Otherwise, no warnings. Check one complete.
Note that there are more people missing on all observations having dropped time points 2,3, and 17. Was 7866, but this version is 8217 .
All resolved. See below sections for descriptions.
Check any errors with running any of the models. Output here will not say which model it belongs to, so will need to manually examine if errors occur,
character(0)
Some of the models have failed to converge:
Steps:
* Increased number of iterations to 5000 (Mplus default is 1000).
* This worked somewhat for the piecewise and quadratic (without pairwise correlations) which now have more informative errors
* This did not work at all for quadratic (with pairwise residual correlation). Increased iterations to 10000 for this one.
* This worked.
Now all three run, but could not compute standard errors as the models were misidentified or did not converge.
Issues are all with specified parameters now:
quadratic: Problem involving parameter 2 (T_3) * examining the data, I have made an error with the dataset. There are 4325 in this time point. This corresponds to the second follow up that is meant to be dropped from this datafile. Interestingly, T4 onwards is correct. T3 should have 2662 observations. Appears I have wrongly dropped this wave instead of follow up 2. Note: this is exactly what I had done. Corrected and rerun all models
Checked this manually, and there are only 4 people who did both timepoint 3 and timepoint 6, 4 who did timepoint 3 and tomepoint 8, and only 3 who did time point 3 and timepoint 16
Options are either dropping time point 3, or dropping coverage requirement.
The two quadratic models are not converging. Without pairwise correlations the parameter MPlus flags is the intercept. With pairwise correlation it flags the correlations between time points 15 and 16
Tabs to check that:
1. Any residual variance is positive. Significant negative residual variance would indicate a model error. 2. The loadings of indicators and intercepts are 1 3. Estimates late correctly to specific time points (i.e. the relative time references for each time point make sense and are accurate). 4. The intercept is zero for all time points 5. Model estimated correlations between time points look sensible
Output from each relevant file will be printed in each tab for review without needing MPlus software. You may need to scroll across to see full output.
check if any residual variances are negative and significant (would indicate model error)
All variances are positive
model <- dep.traj.all$phq_latentGrowthCurve_linear.out$parameters$unstandardized
model[model$paramHeader == "Residual.Variances",]All variances are positive
model <- dep.traj.all$phq_latentGrowthCurve_linear_withpairwiseCorrelations.out$parameters$unstandardized
model[model$paramHeader == "Residual.Variances",]All variances are positive
model <- dep.traj.all$phq_latentGrowthCurve_quadratic.out$parameters$unstandardized
model[model$paramHeader == "Residual.Variances",]All variances are positive
model <- dep.traj.all$phq_latentGrowthCurve_quadratic_withPairwiseCorrelations.out$parameters$unstandardized
model[model$paramHeader == "Residual.Variances",]Check loadings of indicators and intercepts are 1, and estimates relate to specified time points (e.g. lin 0, 1, 2; quad 0, 1, 4)
Note that time references are relative. For most models run here, a difference of 1 refers to a gap of 2 weeks. Thus a difference of 2 refers to a gap of 4 weeks. For Quadratic models, a difference of 0.5 indicates a gap of 2 weeks, and a difference of 1 indicates a gap of 4 weeks. 0 is always the start of a trajectory, before we expect any change to occur. E.g. baseline would be 0 in a linear model.
So linear trajectories (given we are not including follow up time points 1 and 2, which occurred during our 2 week assessment phase) should have a time point reference that looks something like: 0, 3, 4, 5, 6, 7, 8, 9 (two weekly assessment time period ends here, four weeekly assessment period begins), 11, 13, 15, 17, 19, 21, 23, 25
All is fine
dep.traj.all$phq_latentGrowthCurve_linear.out$output[291:326] %>% str_split("#") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
dep.traj.all$phq_latentGrowthCurve_linear_withpairwiseCorrelations.out$output[295:330] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
dep.traj.all$phq_latentGrowthCurve_quadratic.out$output[298:353] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
dep.traj.all$phq_latentGrowthCurve_quadratic_withPairwiseCorrelations.out$output[297:352] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
dep.traj.all$phq_latentGrowthCurve_piecewise.out$output[341:432] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
dep.traj.all$phq_latentGrowthCurve_piecewise_wPairwise_correlations.out$output[345:436] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
check intercept 0 for all time points
All is fine
dep.traj.all$phq_latentGrowthCurve_linear.out$output[427:444] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
dep.traj.all$phq_latentGrowthCurve_linear_withpairwiseCorrelations.out$output[383:400] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
dep.traj.all$phq_latentGrowthCurve_quadratic.out$output[366:383] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
dep.traj.all$phq_latentGrowthCurve_quadratic_withPairwiseCorrelations.out$output[407:423] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
dep.traj.all$phq_latentGrowthCurve_piecewise.out$output[458:475] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
dep.traj.all$phq_latentGrowthCurve_piecewise_wPairwise_correlations.out$output[507:524] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Check that the model estimated correlations between all time points looks sensible.
All is fine
dep.traj.all$phq_latentGrowthCurve_linear.out$output[653:698] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
dep.traj.all$phq_latentGrowthCurve_linear_withpairwiseCorrelations.out$output[746:792] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
dep.traj.all$phq_latentGrowthCurve_quadratic.out$output[709:762] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
dep.traj.all$phq_latentGrowthCurve_quadratic_withPairwiseCorrelations.out$output[792:844] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
dep.traj.all$phq_latentGrowthCurve_piecewise.out$output[854:906] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
dep.traj.all$phq_latentGrowthCurve_piecewise_wPairwise_correlations.out$output[948:1000] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Sorted by BIC, as this is the most informative criteria (as per Muthens)
justSummaries <- do.call("rbind.fill",
sapply(dep.traj.all ,"[", "summaries")) %>%
select(c("Title", "Observations", "Parameters", "AIC", "BIC", "CFI", "TLI", "SRMR", "RMSEA_Estimate")) %>%
arrange(BIC)
justSummaries %>%
kbl() %>%
kable_styling()| Title | Observations | Parameters | AIC | BIC | CFI | TLI | SRMR | RMSEA_Estimate |
|---|---|---|---|---|---|---|---|---|
| Depression symptoms in RAMP and COPING piecewise trajectories With pairwise correlations | 42006 | 46 | 1770170 | 1770568 | 0.995 | 0.994 | 0.009 | 0.015 |
| Depression symptoms in RAMP and COPING Quadratic Growth Curve With pairwise correlations | 42006 | 39 | 1772577 | 1772914 | 0.989 | 0.988 | 0.013 | 0.022 |
| Depression symptoms in RAMP and COPING Linear Growth Curve With pairwise correlations | 42006 | 36 | 1774121 | 1774432 | 0.985 | 0.984 | 0.016 | 0.025 |
| Depression symptoms in RAMP and COPING piecewise trajectory No pairwise correlations | 42006 | 36 | 1774531 | 1774842 | 0.985 | 0.984 | 0.010 | 0.025 |
| Depression symptoms in RAMP and COPING Quadratic Growth Curve No pairwise correlations | 42006 | 25 | 1781344 | 1781560 | 0.968 | 0.970 | 0.017 | 0.035 |
| Depression symptoms in RAMP and COPING Linear Growth Curve No pairwise correlations | 42006 | 21 | 1786100 | 1786282 | 0.956 | 0.960 | 0.021 | 0.040 |
I show the plots for all runs. Output for all runs is available in the GMM folder on GitHub to review input and output data that generated these plots.
Showing model estimated means from piecewise trajectories with pairwise correlations between contiguous time points (solid lines) overlaid on the observed means from our sample (dashed line)
growth_curves_best <- ggplot() +
geom_line(data = dat_dep_best,
aes(x = timepoint, y = score, colour = Model,linetype=Type)) +
geom_point(data = dat_dep_best,
aes(x = timepoint, y = score, colour = Model, shape=Model), show.legend = F) +
scale_color_manual(values = pres.palette) +
labs(x = "Follow up time point", y ="PHQ9 Score") +
scale_y_continuous(expand = c(0,0), limits = c(6,10), breaks=seq(6, 10, 1)) +
scale_x_continuous(expand = c(0,0), limits = c(-0.5,15.5), breaks=seq(0, 15, by = 1),
labels = timepoint_list) +
theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'dashed',
colour = "gray"),
axis.text = element_text(colour="black", size = 14),
axis.title = element_text(colour="black", size = 16),
axis.text.x = element_text(angle = 90,
hjust=1),
legend.key = element_blank(),
legend.text = element_text(colour = "black", size = 11),
panel.background = element_blank())
growth_curves_bestShowing model estimated means from linear, quadratic, and piecewise trajectories (solid lines) overlaid on the observed means from our sample (dashed line) Only showing trajectories that included pairwise correlations between residuals for contiguous time points. This was the best fit for the data.
growth_curves_run2c <- ggplot() +
geom_line(data = dat_dep_run_2c,
aes(x = timepoint, y = score, colour = Model,linetype=Type)) +
geom_point(data = dat_dep_run_2c,
aes(x = timepoint, y = score, colour = Model, shape=Model), show.legend = F) +
scale_color_manual(values = pres.palette) +
labs(x = "Follow up time point", y ="PHQ9 Score") +
scale_y_continuous(expand = c(0,0), limits = c(6,10), breaks=seq(6, 10, 1)) +
scale_x_continuous(expand = c(0,0), limits = c(-0.5,15.5), breaks=seq(0, 15, by = 1),
labels = timepoint_list) +
theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'dashed',
colour = "gray"),
axis.text = element_text(colour="black", size = 14),
axis.title = element_text(colour="black", size = 16),
axis.text.x = element_text(angle = 90,
hjust=1),
legend.key = element_blank(),
legend.text = element_text(colour = "black", size = 11),
panel.background = element_blank())
growth_curves_run2cPlotting the estimated means for the trajectories with adjacent residuals correlated as these all had better fit than uncorrelated models
growth_curves_run2 <- ggplot() +
geom_line(data = dat_dep_run_2,
aes(x = timepoint, y = score, colour = Model), size = 1) +
geom_point(data = dat_dep_run_2,
aes(x = timepoint, y = score, colour = Model, shape=Model), show.legend = F) +
scale_color_manual(values = pres.palette) +
labs(x = "Follow up time point", y ="PHQ9 Score") +
scale_y_continuous(expand = c(0,0), limits = c(6,10), breaks=seq(6, 10, 1)) +
scale_x_continuous(expand = c(0,0), limits = c(-0.5,15.5), breaks=seq(0, 15, by = 1),
labels = timepoint_list) +
theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'dashed',
colour = "gray"),
axis.text = element_text(colour="black", size = 14),
axis.title = element_text(colour="black", size = 16),
axis.text.x = element_text(angle = 90,
hjust=1),
legend.key = element_blank(),
legend.text = element_text(colour = "black", size = 11),
panel.background = element_blank())
growth_curves_run2Show model that performs best of each pair of models (with or without residual correlations between contiguous data points)
Linear: with correlations quadratic: with correlations Piecewise: with correlations
growth_curves_run2b <- ggplot() +
geom_line(data = dat_dep_run_2b,
aes(x = timepoint, y = score, colour = Model), size = 1) +
geom_point(data = dat_dep_run_2b,
aes(x = timepoint, y = score, colour = Model, shape=Model), show.legend = F) +
scale_color_manual(values = pres.palette[c(4, 1:3)]) +
labs(x = "Follow up time point", y ="PHQ9 Score") +
scale_y_continuous(expand = c(0,0), limits = c(6,10), breaks=seq(6, 10, 1)) +
scale_x_continuous(expand = c(0,0), limits = c(-0.5,15.5), breaks=seq(0, 15, by = 1),
labels = timepoint_list) +
theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'dashed',
colour = "gray"),
axis.text = element_text(colour="black", size = 14),
axis.title = element_text(colour="black", size = 16),
axis.text.x = element_text(angle = 90,
hjust=1),
legend.key = element_blank(),
legend.text = element_text(colour = "black", size = 11),
panel.background = element_blank())
growth_curves_run2bgad_check <- readRDS("../../../data_clean/gad/gad.clean_merged_total_scores.rds")
# get a list of the total score columns
total_cols <- grep("total_Wave", names(gad_check))
#count how many columns we have
n_cols <- length(total_cols)
# create a list of how many NAs are missing across these 18 columns per row
list_n_missing <- apply(gad_check[total_cols],1,function(x) sum(is.na(x)))
# count how many people are missing data for all columns (NA = total_cols)
list_all_missing <- sum(list_n_missing ==n_cols)
print(dim(gad_check)[1])[1] 52414
[1] 10968
this reads in all output, including the gh5 plot files
at the moment this doesnt say which model the warning is from, so any warning requires checking to identify its source. Though order in list will give some indication
[1] “c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10969", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10969", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("All continuous latent variable covariances involving LIN4 have been fixed to 0", "because the variance of LIN4 is fixed at 0.") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10969", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10969", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10969", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS")”
[2] “c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10969", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10969", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10969", "2 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10969", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10969", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10969", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS")”
Warnings are all about missing data (missing on all obs, this number lines up with our expectations from the dataset in the section above)
Otherwise, no warnings. Check one complete.
No errors occurred. Used the same alterations as for PHQ (see above)
Check any errors with running any of the models. Output here will not say which model it belongs to, so will need to manually examine if errors occur,
character(0)
Tabs to check that:
1. Any residual variance is positive. Significant negative residual variance would indicate a model error. 2. The loadings of indicators and intercepts are 1 3. Estimates late correctly to specific time points (i.e. the relative time references for each time point make sense and are accurate). 4. The intercept is zero for all time points 5. Model estimated correlations between time points look sensible
Output from each relevant file will be printed in each tab for review without needing MPlus software. You may need to scroll across to see full output.
check if any residual variances are negative and significant (would indicate model error) ###### Linear All variances are positive
model <- anx.traj.all$gad_latentGrowthCurve_linear.out$parameters$unstandardized
model[model$paramHeader == "Residual.Variances",]All variances are positive
model <- anx.traj.all$gad_latentGrowthCurve_linear_withpairwiseCorrelations.out$parameters$unstandardized
model[model$paramHeader == "Residual.Variances",]All variances are positive
model <- anx.traj.all$gad_latentGrowthCurve_quadratic.out$parameters$unstandardized
model[model$paramHeader == "Residual.Variances",]All variances are positive
model <- anx.traj.all$gad_latentGrowthCurve_quadratic_withPairwiseCorrelations.out$parameters$unstandardized
model[model$paramHeader == "Residual.Variances",]Check loadings of indicators and intercepts are 1, and estimates relate to specified time points (e.g. lin 0, 1, 2; quad 0, 1, 4)
Note that time references are relative. For most models run here, a difference of 1 refers to a gap of 2 weeks. Thus a difference of 2 refers to a gap of 4 weeks. For Quadratic models, a difference of 0.5 indicates a gap of 2 weeks, and a difference of 1 indicates a gap of 4 weeks. 0 is always the start of a trajectory, before we expect any change to occur. E.g. baseline would be 0 in a linear model.
So linear trajectories (given we are not including follow up time points 1 and 2, which occurred during our 2 week assessment phase) should have a time point reference that looks something like: 0, 3, 4, 5, 6, 7, 8, 9 (two weekly assessment time period ends here, four weeekly assessment period begins), 11, 13, 15, 17, 19, 21, 23, 25
All is fine
anx.traj.all$gad_latentGrowthCurve_linear.out$output[290:327] %>%
str_split("/t") %>%
kbl() %>%
kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anx.traj.all$gad_latentGrowthCurve_linear_withpairwiseCorrelations.out$output[293:330] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anx.traj.all$gad_latentGrowthCurve_quadratic.out$output[298:353] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anx.traj.all$gad_latentGrowthCurve_quadratic_withPairwiseCorrelations.out$output[297:352] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anx.traj.all$gad_latentGrowthCurve_piecewise.out$output[340:432] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anx.traj.all$gad_latentGrowthCurve_piecewise_wPairwise_correlations.out$output[558:648] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
check intercept 0 for all time points
All is fine
anx.traj.all$gad_latentGrowthCurve_linear.out$output[335:352] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anx.traj.all$gad_latentGrowthCurve_linear_withpairwiseCorrelations.out$output[383:400] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anx.traj.all$gad_latentGrowthCurve_quadratic.out$output[366:383] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anx.traj.all$gad_latentGrowthCurve_quadratic_withPairwiseCorrelations.out$output[410:427] %>%
str_split("/t") %>%
kbl() %>%
kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anx.traj.all$gad_latentGrowthCurve_piecewise.out$output[458:475] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anx.traj.all$gad_latentGrowthCurve_piecewise_wPairwise_correlations.out$output[718:736] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Check that the model estimated correlations between all time points looks sensible.
All is fine
anx.traj.all$gad_latentGrowthCurve_linear.out$output[653:705] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anx.traj.all$gad_latentGrowthCurve_linear_withpairwiseCorrelations.out$output[746:798] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anx.traj.all$gad_latentGrowthCurve_quadratic.out$output[709:760] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anx.traj.all$gad_latentGrowthCurve_quadratic_withPairwiseCorrelations.out$output[851:798] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anx.traj.all$gad_latentGrowthCurve_piecewise.out$output[854:906] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anx.traj.all$gad_latentGrowthCurve_piecewise_wPairwise_correlations.out$output[948:1000] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Sorted by BIC, as this is the most informative criteria (as per Muthens)
justSummaries <- do.call("rbind.fill",
sapply(anx.traj.all ,"[", "summaries")) %>%
select(c("Title", "Observations", "Parameters", "AIC", "BIC", "CFI", "TLI", "SRMR", "RMSEA_Estimate")) %>%
arrange(BIC)
justSummaries %>%
kbl() %>%
kable_styling()| Title | Observations | Parameters | AIC | BIC | CFI | TLI | SRMR | RMSEA_Estimate |
|---|---|---|---|---|---|---|---|---|
| Anxiety symptoms in RAMP and COPING piecewise trajectory With pairwise correlations | 41445 | 46 | 1658013 | 1658410 | 0.996 | 0.996 | 0.008 | 0.013 |
| Anxiety symptoms in RAMP and COPING Quadratic Growth Curve With pairwise correlations | 41445 | 40 | 1660647 | 1660992 | 0.989 | 0.988 | 0.013 | 0.021 |
| Anxiety symptoms in RAMP and COPING piecewise trajectory No pairwise correlations | 41445 | 36 | 1661334 | 1661645 | 0.988 | 0.987 | 0.009 | 0.022 |
| Anxiety symptoms in RAMP and COPING Linear Growth Curve With pairwise correlations | 41445 | 36 | 1661947 | 1662258 | 0.986 | 0.985 | 0.016 | 0.024 |
| Anxiety symptoms in RAMP and COPING Quadratic Growth Curve No pairwise correlations | 41445 | 25 | 1667841 | 1668056 | 0.971 | 0.972 | 0.018 | 0.032 |
| Anxiety symptoms in RAMP and COPING Linear Growth Curve No pairwise correlations | 41445 | 21 | 1671934 | 1672116 | 0.960 | 0.963 | 0.021 | 0.037 |
I show the plots for all runs. Output for all runs is available in the GMM folder on GitHub to review input and output data that generated these plots.
Showing model estimated means from piecewise trajectories with pairwise correlations between contiguous time points (solid lines) overlaid on the observed means from our sample (dashed line)
growth_curves_best <- ggplot() +
geom_line(data = dat_anx_best,
aes(x = timepoint, y = score, colour = Model,linetype=Type)) +
geom_point(data = dat_anx_best,
aes(x = timepoint, y = score, colour = Model, shape=Model), show.legend = F) +
scale_color_manual(values = pres.palette) +
labs(x = "Follow up time point", y ="GAD7 Score") +
scale_y_continuous(expand = c(0,0), limits = c(5,7), breaks=seq(5, 7, 1)) +
scale_x_continuous(expand = c(0,0), limits = c(-0.5,15.5), breaks=seq(0, 15, by = 1),
labels = timepoint_list) +
theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'dashed',
colour = "gray"),
axis.text = element_text(colour="black", size = 14),
axis.title = element_text(colour="black", size = 16),
axis.text.x = element_text(angle = 90,
hjust=1),
legend.key = element_blank(),
legend.text = element_text(colour = "black", size = 11),
panel.background = element_blank())
growth_curves_best #### Observed and estimated means
Plot estimated means from model (solid lines) alongside the observed means (dashed lines) from the data.
Using the best models from each model type pairing (i.e. the version with pairwise correlations of residuals between contiguous time points)
growth_curves_run2c <- ggplot() +
geom_line(data = dat_anx_run_2c,
aes(x = timepoint, y = score, colour = Model,linetype=Type)) +
geom_point(data = dat_anx_run_2c,
aes(x = timepoint, y = score, colour = Model, shape=Model), show.legend = F) +
scale_color_manual(values = pres.palette) +
labs(x = "Follow up time point", y ="GAD7 Score") +
scale_y_continuous(expand = c(0,0), limits = c(5,7), breaks=seq(5, 7, 1)) +
scale_x_continuous(expand = c(0,0), limits = c(-0.5,15.5), breaks=seq(0, 15, by = 1),
labels = timepoint_list) +
theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'dashed',
colour = "gray"),
axis.text = element_text(colour="black", size = 14),
axis.title = element_text(colour="black", size = 16),
axis.text.x = element_text(angle = 90,
hjust=1),
legend.key = element_blank(),
legend.text = element_text(colour = "black", size = 11),
panel.background = element_blank())
growth_curves_run2cPlotting the estimated means for the trajectories with adjacent residuals correlated as these all had better fit than uncorrelated models
growth_curves_run2 <- ggplot() +
geom_line(data = dat_anx_run_2,
aes(x = timepoint, y = score, colour = Model), size = 1) +
geom_point(data = dat_anx_run_2,
aes(x = timepoint, y = score, colour = Model, shape=Model), show.legend = F) +
scale_color_manual(values = pres.palette) +
labs(x = "Follow up time point", y ="GAD7 Score") +
scale_y_continuous(expand = c(0,0), limits = c(5,7), breaks=seq(5, 7, 1)) +
scale_x_continuous(expand = c(0,0), limits = c(-0.5,15.5), breaks=seq(0, 15, by = 1),
labels = timepoint_list) +
theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'dashed',
colour = "gray"),
axis.text = element_text(colour="black", size = 14),
axis.title = element_text(colour="black", size = 16),
axis.text.x = element_text(angle = 90,
hjust=1),
legend.key = element_blank(),
legend.text = element_text(colour = "black", size = 11),
panel.background = element_blank())
growth_curves_run2Only show the best of each pair of models (with or without pairwise residual correlations)
Linear: With correlations Quadratic: With correlations Piecewise: With correlations
growth_curves_run2b <- ggplot() +
geom_line(data = dat_anx_run_2b,
aes(x = timepoint, y = score, colour = Model), size = 1) +
geom_point(data = dat_anx_run_2b,
aes(x = timepoint, y = score, colour = Model, shape=Model), show.legend = F) +
scale_color_manual(values = pres.palette[c(4, 1:3)]) +
labs(x = "Follow up time point", y ="GAD7 Score") +
scale_y_continuous(expand = c(0,0), limits = c(5,7), breaks=seq(5, 7, 1)) +
scale_x_continuous(expand = c(0,0), limits = c(-0.5,15.5), breaks=seq(0, 15, by = 1),
labels = timepoint_list) +
theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'dashed',
colour = "gray"),
axis.text = element_text(colour="black", size = 14),
axis.title = element_text(colour="black", size = 16),
axis.text.x = element_text(angle = 90,
hjust=1),
legend.key = element_blank(),
legend.text = element_text(colour = "black", size = 11),
panel.background = element_blank())
growth_curves_run2bMeasure has been reverse coded, such that higher scores indicate higher levels of negative affect, to correspond with the direction of the GAD-7 and PHQ 9
Will run these trajectories starting from follow up time 3 as this is the first point where we have both CPOING and RAMP data
masq_check <- readRDS("../../../data_clean/masq/masq.clean_merged_total_scores.rds")
# get a list of the total score columns
total_cols <- grep("total_Wave", names(masq_check))
#count how many columns we have
n_cols <- length(total_cols)
# create a list of how many NAs are missing across these 18 columns per row
list_n_missing <- apply(masq_check[total_cols],1,function(x) sum(is.na(x)))
# count how many people are missing data for all columns (NA = total_cols)
list_all_missing <- sum(list_n_missing ==n_cols)
print(dim(masq_check)[1])[1] 37151
[1] 1320
Note: this corresponds to the dataset with all time points, thus, different missing on all data pattern from MPlus run with time points 1,2 and baseline (have more missing all when we exclude the RAMP baseline!)
this reads in all output, including the gh5 plot files
at the moment this doesnt say which model the warning is from, so any warning requires checking to identify its source. Though order in list will give some indication
[1] “c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("All continuous latent variable covariances involving LIN1 have been fixed to 0", "because the variance of LIN1 is fixed at 0.") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS")”
[2] “c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("All continuous latent variable covariances involving LIN4 have been fixed to 0", "because the variance of LIN4 is fixed at 0.") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS")”
[3] “c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "3 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS")”
Warnings are all about missing data (missing on all obs, this number lines up with our expectations from the dataset in the section above)
Otherwise, no warnings. Check one complete.
Check any errors with running any of the models. Output here will not say which model it belongs to, so will need to manually examine if errors occur,
character(0)
piecewise with correlations did not run due to issues with the first linear piece. Likely due to few pieces contributing to this section. Restrained variance to 0 for this piece and it ran successfully.
Tabs to check that:
1. Any residual variance is positive. Significant negative residual variance would indicate a model error. 2. The loadings of indicators and intercepts are 1 3. Estimates late correctly to specific time points (i.e. the relative time references for each time point make sense and are accurate). 4. The intercept is zero for all time points 5. Model estimated correlations between time points look sensible
Output from each relevant file will be printed in each tab for review without needing MPlus software. You may need to scroll across to see full output.
check if any residual variances are negative and significant (would indicate model error) ###### Linear All variances are positive
model <- anhed.traj.all$masq_latentGrowthCurve_linear.out$parameters$unstandardized
model[model$paramHeader == "Residual.Variances",]All variances are positive
model <- anhed.traj.all$masq_latentGrowthCurve_linear_withpairwiseCorrelations.out$parameters$unstandardized
model[model$paramHeader == "Residual.Variances",]All variances are positive
model <- anhed.traj.all$masq_latentGrowthCurve_quadratic.out$parameters$unstandardized
model[model$paramHeader == "Residual.Variances",]All variances are positive
model <- anhed.traj.all$masq_latentGrowthCurve_quadratic_withPairwiseCorrelations.out$parameters$unstandardized
model[model$paramHeader == "Residual.Variances",]Check loadings of indicators and intercepts are 1, and estimates relate to specified time points (e.g. lin 0, 1, 2; quad 0, 1, 4)
Note that time references are relative. For most models run here, a difference of 1 refers to a gap of 2 weeks. Thus a difference of 2 refers to a gap of 4 weeks. For Quadratic models, a difference of 0.5 indicates a gap of 2 weeks, and a difference of 1 indicates a gap of 4 weeks. 0 is always the start of a trajectory, before we expect any change to occur. E.g. baseline would be 0 in a linear model.
So linear trajectories (given we are not including follow up time points 1 and 2, which occurred during our 2 week assessment phase) should have a time point reference that looks something like: 0, 3, 4, 5, 6, 7, 8, 9 (two weekly assessment time period ends here, four weeekly assessment period begins), 11, 13, 15, 17, 19, 21, 23, 25
All is fine
anhed.traj.all$masq_latentGrowthCurve_linear.out$output[280:314] %>%
str_split("/t") %>%
kbl() %>%
kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anhed.traj.all$masq_latentGrowthCurve_linear_withpairwiseCorrelations.out$output[284:317] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anhed.traj.all$masq_latentGrowthCurve_quadratic.out$output[400:450] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anhed.traj.all$masq_latentGrowthCurve_quadratic_withPairwiseCorrelations.out$output[287:337] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anhed.traj.all$masq_latentGrowthCurve_piecewise.out$output[490:574] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anhed.traj.all$masq_latentGrowthCurve_piecewise_wPairwise_correlations.out$output[340:424] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
check intercept 0 for all time points
All is fine
anhed.traj.all$masq_latentGrowthCurve_linear.out$output[322:338] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anhed.traj.all$masq_latentGrowthCurve_linear_withpairwiseCorrelations.out$output[367:383] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anhed.traj.all$masq_latentGrowthCurve_quadratic.out$output[463:479] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anhed.traj.all$masq_latentGrowthCurve_quadratic_withPairwiseCorrelations.out$output[392:408] %>%
str_split("/t") %>%
kbl() %>%
kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anhed.traj.all$masq_latentGrowthCurve_piecewise.out$output[600:616] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anhed.traj.all$masq_latentGrowthCurve_piecewise_wPairwise_correlations.out$output[694:710] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Check that the model estimated correlations between all time points looks sensible.
All is fine
anhed.traj.all$masq_latentGrowthCurve_linear.out$output[600:643] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anhed.traj.all$masq_latentGrowthCurve_linear_withpairwiseCorrelations.out$output[687:730] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anhed.traj.all$masq_latentGrowthCurve_quadratic.out$output[654:697] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anhed.traj.all$masq_latentGrowthCurve_quadratic_withPairwiseCorrelations.out$output[736:779] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anhed.traj.all$masq_latentGrowthCurve_piecewise.out$output[793:836] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
All is fine
anhed.traj.all$masq_latentGrowthCurve_piecewise_wPairwise_correlations.out$output[887:930] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Sorted by BIC, as this is the most informative criteria (as per Muthens)
justSummaries <- do.call("rbind.fill",
sapply(anhed.traj.all ,"[", "summaries")) %>%
select(c("Title", "Observations", "Parameters", "AIC", "BIC", "CFI", "TLI", "SRMR", "RMSEA_Estimate")) %>%
arrange(BIC)
justSummaries %>%
kbl() %>%
kable_styling()| Title | Observations | Parameters | AIC | BIC | CFI | TLI | SRMR | RMSEA_Estimate |
|---|---|---|---|---|---|---|---|---|
| Anhedonia symptoms in RAMP and COPING piecewise trajectory With pairwise correlations | 32892 | 40 | 1795036 | 1795372 | 0.984 | 0.982 | 0.032 | 0.032 |
| Anhedonia symptoms in RAMP and COPING piecewise trajectory No pairwise correlations | 32892 | 35 | 1796887 | 1797181 | 0.978 | 0.977 | 0.026 | 0.037 |
| Anhedonia symptoms in RAMP and COPING Quadratic Growth Curve With pairwise correlations | 32892 | 38 | 1797427 | 1797746 | 0.975 | 0.973 | 0.042 | 0.040 |
| Anhedonia symptoms in RAMP and COPING Linear Growth Curve With pairwise correlations | 32892 | 34 | 1799404 | 1799690 | 0.968 | 0.967 | 0.043 | 0.044 |
| Anhedonia symptoms in RAMP and COPING Quadratic Growth Curve No pairwise correlations | 32892 | 24 | 1802466 | 1802668 | 0.958 | 0.961 | 0.044 | 0.048 |
| Anhedonia symptoms in RAMP and COPING Linear Growth Curve No pairwise correlations | 32892 | 20 | 1806998 | 1807166 | 0.943 | 0.948 | 0.046 | 0.056 |
I show the plots for all runs. Output for all runs is available in the GMM folder on GitHub to review input and output data that generated these plots.
Showing model estimated means from piecewise trajectories with pairwise correlations between contiguous time points (solid lines) overlaid on the observed means from our sample (dashed line)
growth_curves_best <- ggplot() +
geom_line(data = dat_anhed_best,
aes(x = timepoint, y = score, colour = Model,linetype=Type)) +
geom_point(data = dat_anhed_best,
aes(x = timepoint, y = score, colour = Model, shape=Model), show.legend = F) +
scale_color_manual(values = pres.palette) +
labs(x = "Follow up time point", y ="masq (reverse) Score") +
scale_y_continuous(expand = c(0,0), limits = c(25,29), breaks=seq(25, 29, 1)) +
scale_x_continuous(expand = c(0,0), limits = c(-0.5,14.5), breaks=seq(0, 14, by = 1),
labels = timepoint_list_anh) +
theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'dashed',
colour = "gray"),
axis.text = element_text(colour="black", size = 14),
axis.title = element_text(colour="black", size = 16),
axis.text.x = element_text(angle = 90,
hjust=1),
legend.key = element_blank(),
legend.text = element_text(colour = "black", size = 11),
panel.background = element_blank())
growth_curves_best #### Observed and estimated means
Plot estimated means from model (solid lines) alongside the observed means (dashed lines) from the data.
Using the best models from each model type pairing (i.e. the version with pairwise correlations of residuals between contiguous time points)
growth_curves_run2c <- ggplot() +
geom_line(data = dat_anhed_run_2c,
aes(x = timepoint, y = score, colour = Model,linetype=Type)) +
geom_point(data = dat_anhed_run_2c,
aes(x = timepoint, y = score, colour = Model, shape=Model), show.legend = F) +
scale_color_manual(values = pres.palette) +
labs(x = "Follow up time point", y ="masq (reverse) Score") +
scale_y_continuous(expand = c(0,0), limits = c(25,29), breaks=seq(25, 29, 1)) +
scale_x_continuous(expand = c(0,0), limits = c(-0.5,14.5), breaks=seq(0, 14, by = 1),
labels = timepoint_list_anh) +
theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'dashed',
colour = "gray"),
axis.text = element_text(colour="black", size = 14),
axis.title = element_text(colour="black", size = 16),
axis.text.x = element_text(angle = 90,
hjust=1),
legend.key = element_blank(),
legend.text = element_text(colour = "black", size = 11),
panel.background = element_blank())
growth_curves_run2cPlotting the estimated means for the trajectories with adjacent residuals correlated as these all had better fit than uncorrelated models
growth_curves_run2 <- ggplot() +
geom_line(data = dat_anhed_run_2,
aes(x = timepoint, y = score, colour = Model), size = 1) +
geom_point(data = dat_anhed_run_2,
aes(x = timepoint, y = score, colour = Model, shape=Model), show.legend = F) +
scale_color_manual(values = pres.palette) +
labs(x = "Follow up time point", y ="masq (reverse) Score") +
scale_y_continuous(expand = c(0,0), limits = c(25,29), breaks=seq(25, 29, 1)) +
scale_x_continuous(expand = c(0,0), limits = c(-0.5,14.5), breaks=seq(0, 14, by = 1),
labels = timepoint_list_anh) +
theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'dashed',
colour = "gray"),
axis.text = element_text(colour="black", size = 14),
axis.title = element_text(colour="black", size = 16),
axis.text.x = element_text(angle = 90,
hjust=1),
legend.key = element_blank(),
legend.text = element_text(colour = "black", size = 11),
panel.background = element_blank())
growth_curves_run2Only show the best of each pair of models (with or without pairwise residual correlations)
Linear: With correlations Quadratic: With correlations Piecewise: With correlations
growth_curves_run2b <- ggplot() +
geom_line(data = dat_anhed_run_2b,
aes(x = timepoint, y = score, colour = Model), size = 1) +
geom_point(data = dat_anhed_run_2b,
aes(x = timepoint, y = score, colour = Model, shape=Model), show.legend = F) +
scale_color_manual(values = pres.palette[c(4, 1:3)]) +
labs(x = "Follow up time point", y ="masq (reverse) Score") +
scale_y_continuous(expand = c(0,0), limits = c(25,29), breaks=seq(25, 29, 1)) +
scale_x_continuous(expand = c(0,0), limits = c(-0.5,14.5), breaks=seq(0, 14, by = 1),
labels = timepoint_list_anh) +
theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'dashed',
colour = "gray"),
axis.text = element_text(colour="black", size = 14),
axis.title = element_text(colour="black", size = 16),
axis.text.x = element_text(angle = 90,
hjust=1),
legend.key = element_blank(),
legend.text = element_text(colour = "black", size = 11),
panel.background = element_blank())
growth_curves_run2b